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Abstract. We propose a Markov model for the spread of Hepatitis C virus 
(HCV) among drug users who use injections. We then proceed to an asymp- 
totic analysis (large initial population) and show that the Markov process is 
close to the solution of a non linear autonomous differential system. We prove 
both a law of large numbers and functional central limit theorem to precise 
the speed of convergence towards the limiting system. The deterministic sys- 
tem itself converges, as time goes to infinity, to an equilibrium point. This 
corroborates the empirical observations about the prevalence of HCV. 



1. Motivations 

Hepatitis C virus (HCV) infects 170 million people in the world (3 % of the 
population) and 9 million in Europe (1 % of the population) [16]. More than 75 % 
of newly infected patients progress to develop chronic infection. Then, Cirrhosis 
develops in about 10 % to 20 %, and liver cancer develops in 1 % to 5 % over a 
period of 20 to 30 years. These long-term consequences, which suggest an increased 
mortality due to HCV infection, make the prevention of spread of hepatitis C a 
major public health concern. 

HCV is spread primarily by direct contact with human blood. In developed 
countries that have safe blood supplies, the population infected by HCV is closely 
related to injecting drug users (IDU). It is estimated that 90 % of infectious are due 
to IDU [11]. In order to reduce the numbers of new hepatitis C cases, preventing 
infections in IDU is then a priority. Programs exist all over the world which try to 
reduce the prevalence of many infectious diseases like HIV or hepatitis C, among 
injecting drug users. They are mainly based on needle exchanges. It turns out that 
after several years of such programs, the HIV prevalence seems to be now rather 
low whereas the percentage of IDU who are HCV positive remains about 60 % 
[9, 11]. We were asked by epidemiologists to provide them a mathematical model 
which could quantitatively evaluate the differences between the two diseases. 

It is always a challenge to analyze an epidemic problem because there are so 
many real-life situations that should be incorporated while keeping the mathemat- 
ical model tractable. Moreover, epidemic field studies are expensive and hard to 
organize so that parameter estimates are rare and often imprecise. It is thus nec- 
essary to deal with parsimonious models whose parameters have clear and visible 
meaning. To the best of our knowledge, the only models which have been developed 
for the dynamics of HCV transmission are found in the references [14, 4] . It is a de- 
terministic model with more than twenty-five parameters, for which the authors do 
not have explicit results for the asymptotics and only estimate them by simulations. 
In our paper, we propose a parsimonious Markovian model for the spread of HCV 
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in a local population of IDU. It should be noted that our model bears some resem- 
blance to a random SIR (Susceptible-Infected-Recovered) model but differs from it 
by some essential characteristics. Our Susceptible (respectively Infected) are IDU 
who are sero- negative (respectively sero-positive) . There is no Recovered category 
in our model since we can't measure their number (when they are no longer IDU, 
they can't be counted in studies focused on drugs users). Moreover, our population 
is not closed (there are new susceptible all the time) and a new drug user may be 
infected at his first injection. This means that there is an exogeneous flow to the 
Infected category, a feature which is not included in usual SIR models. 

To keep the Markovian character of our model, we made the following usual 
and reasonable hypothesis. Exogenous antibody-positive and antibody-negative 
individuals arrive in this local population according to Poisson processes. If initiated 
by an antibody-positive drug addict, a new IDU acquires the virus very rapidly after 
the initiation [1, 6]. HCV then spreads in the population by sharing syringe, needles 
and other accessories (cotton, boilers, etc.). Each individual of the population 
stays in his state (infected/non infected) for an exponentially distributed time. We 
present the model in Section 3. If we denote by Xi(t) (resp. X 2 (t)) the number of 
antibody-positive (resp. antibody-negative) individuals in the local population at 
time t, we prove that the process X = (Xi,X 2 ) is an ergodic Markov process. In 
Section 4., we give a related deterministic differential system connected with this 
Markov process. We study its asymptotic behaviour and give an explicit expression 
of the limit of the solution. In Section 5., we give a mean-field approximation of 
the process X: For large populations, we prove that the process X is close to the 
solution ip of the deterministic differential system. In Section 6., we prove that, 
for large populations, the invariant distribution for the Markov process X can be 
approximated by the Dirac measure which only charges ip(oo). Hence we can give 
an explicit limit of the prevalence of HCV in the population. In Section 7, we 
give a central limit theorem for the approximation of X by ip when the population 
tends to infinity. In Section 8, we show that even for a small value of N, there is a 
good accordance between the prevalence computed on the deterministic limit and 
the prevalence observed on the stochastic model. We also show that this can be 
extended to the sensitivity of the model with respect to slight variations of some 
parameters. 

2. Preliminaries 

Let us denote by ©([0, T], R 2 ) the set of cadlag functions equipped with its usual 
topology. In this Section, we recall some results about cadlag semi-martingales; for 
details we refer to [10]. We assume that we are given (Q, (F tl t > 0), P) a filtered 
probability space satisfying the so-called usual hypothesis. On (fi, (.F t , t> 0), P), 
let X and Y be two real-valued cadlag square integrable semi-martingales. The 
mutual variation of X and Y, denoted by [X, Y], is the right continuous process 
with finite variation such that the following integration by parts formula is satisfied: 

X(t)Y(t)-X(0)Y(0)= [ X(s-)dY(s)+ [ Y(s-)dX(s) + [X,Y] t . 
J(o,t] J(o,t] 

The Meyer process of the couple (X, Y), or its square bracket, is denoted by ( X, Y) 
is the unique right continuous with finite variation predictable process such that 

X(t)Y(t)-X(0)Y(0)-(X, Y) t 

is a martingale. Alternatively, (X,Y) and is the unique right continuous, pre- 
dictable with finite variation, process such that [X, Y] — (X, Y) is a martingale. 
For a vector valued semi- martingale X = {X\, X 2 ) where X\ and X 2 are real valued 
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martingales, we denote by ({X}}, its square bracket, denned by 




In the sequel, if a; is a vector (resp. M a matrix) we denote by ||x|| (resp. ||M||) its 
L 1 -norm. 

Let E be a discrete denumerable space. Let (X(t), t > 0) be an E- valued, pure 
jump Markov process, with infinitesimal generator Q = (q xy , (x, y) £ E x E). For 
any F : B^R, Dynkin's Lemma states that the process: 



Here and hereafter, we identify the matrix Q and the operator Q defined as above. 



We consider the dynamics of HCV among a local population which suffers a con- 
tinuous arrival of exogenous antibody-positive individuals, described by a Poisson 
process of intensity r. We let X\(t) and X2(t) denote the number of antibody- 
positive, respectively antibody-negative, users at time t in the population under 
consideration. The new susceptible drug users arrive as a Poisson process of inten- 
sity A. We assume that for their first injection, they arc initiated by an older IDU 
who has a probability q(t) = Xi(t){Xi{t)+X 2 (t))^ 1 of being infected. For different 
reasons, even in this situation, the probability of being infected, is not exactly one 
and is denoted by pj. Each time, an antibody-negative IDU has an injection, he 
may share some of his paraphernalia and may become infected if the sharing occurs 
with an infected IDU. We summarize all these probabilities by saying that at each 
injection, the probability of becoming infected is pq(t), where p is a parameter to 
be estimated, as is pj. If we denote by a the rate at which an IDU injects, and if 
ap is small, we can assume that the rate at which a sane IDU in the population is 
infected, is given by apq(t). Once infected, an IDU may exit from the population 
under consideration either by a death, self healing or stopping drug usage. The 
whole of these situations is modeled by an exponentially distributed duration with 
parameter p,\. For antibody- negative IDU, the only way to exit the population is 
by stopping drug injection, supposed to happen after an exponentially distributed 
duration with parameter /i 2 • In summary, the transitions are described in Figure 1 . 

For further references, we set 




is a local martingale, where 



QF(x) = ^(F(y)-F(x))q.. 



3. Markov model 



qi(m, n 2 ) = r + Xpi 



m + n 2 



q 2 (n Xl n 2 ) = Mi n\ 
93 n 2 ) = apn 2 



m + n 2 
m 



94 (m, n 2 ) = A(l - pi 
95 (ni, n 2 ) = fi 2 n 2 . 



m + n 2 



) 



Lemma 3.1. Let x° = (x^x®)- Conditionally on X(0) = x° , the process W(t) = 
Xi(t) + X 2 (t) — (xi +x 2 ) is dominated (for the strong stochastic order of processes) 
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ni 

apri2 

ni + n 2 



A(l-pj 




Figure 1 . Transitions of the Markov model. 



6y a Poisson process of intensity r + A. In particular, for any t € [0, T], 



E 



sup||X(i 


)ir 


X(0) = .T° 


t<T 







<(||x°|| + (r + A)rf, 



/or any p > 1 . 



Proof. It suffices to say that by suppressing all the departures, we get another 
system with a population larger than that of the system under consideration, at 
any time, for any trajectory. Then, X x (t) + X 2 if) — (x\ + x^) is less than the 
number of arrivals of a Poisson process of intensity r + X. Since a Poisson process 
has increasing path, its supremum over [0, T] is its value at time T. The second 
assertion follows. □ 

Theorem 3.1. The Markov process X = (Xi, X2) is ergodic. For r > 0, the 
process X is irreducible. For r — 0, the set {(ni, n-i) G N x N, m = 0} is a proper 
closed subset. 

Proof. Let S be the function defined on N x N by 

S(m, n 2 ) = ||(ni,n 2 )|| = m +n 2 . 
If we denote by Q the infinitesimal generator of X, we have 

QS(ni, n 2 ) = X + r - - ^ 2 n 2 - 

Let if be a real strictly greater than (A + r + I)/^_ where /x_ = fiiAfi 2 and consider 
the following finite subset of the state space: 

F> K = {(ni, n 2 ) e N x N, n x + n 2 < K}. 

If (m, n 2 ) belongs to D C K , then 

QS(ni, n 2 ) < X + r - + n 2 ) < -I. 
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Lemma 3.1 implies that both 



E 



sup S(X(s)) 

56 [0, 1] 



and E 



f 1 \QS(X(s))\ds 
Jo 



are finite. Then according to [12, Proposition 8.14], X is ergodic. 

The second and third assertions are immediate through inspection of the transi- 
tion rates. □ 

With the non-linearity appearing in the transitions, it seems hopeless to find 
an exact expression for the stationary probability of the Markov process (Xi, X 2 ). 
As usual in queueing theory [12], we then resort to asymptotic analysis in order 
to gain some insights on the evolution of this system. This means that we let the 
initial population becoming larger and larger. For keeping other quantities of the 
same order of magnitude, one are thus led to increase r and A at the same speed, 
i.e., keeping the ratio i = (r + X)/(xi + xi) constant. Note that in epidemiological 
language, i is the incidence of new susceptible. It is measured in percentage of 
individuals per unit of year. 

4. A DETERMINISTIC DIFFERENTIAL SYSTEM 

The mean field approximation will lead us to investigate the solutions of the fol- 
lowing differential system with initial condition a; = (x^x®) € (R + xR + )\ {(0, 0)}: 



(Sr(x°)) 



V'i(o) 



= r + Xpi 



= A(l-p/ 



- fix ip!(t) + ap 



) - M2 ip2{t) - ap 



Theorem 4.1. For any x° = {x\,x%) € (R+ x R+)\ {(0, 0)} ; there exists a unique 
solution to (S r (x )). Furthermore, this solution is defined on R. For r > 0, the 
differential system has a unique fixed point (£i, £2) in R+ x R-h defined by the 
equations 

(1) £2 = —{r + A - and £1 = — 



c + sgn(a)y/(ab — c) 2 + Aabrfii 



p 2 2ajUi 
where a = ap — Hi + ^2, b = r + X and c = r\i\ + A(l — pi)p 2 - Moreover, for r > 
and any x" G R+\{0, 0}, 



lim (Mt),Mt)) = (Zi, 



Ifr = and x\ = then 



^(t) = for all t and lim (ipi(t),ip2{t)) = (0, A//U 2 ). 

t— 7- + 0O 

If r = and p = ap + ji 2 pi — Ml > 0, then there exists two equilibrium points: 
one is (0, A//U2) and the other is the unique solution with positive first coordinate 
of (1). Ifx\>Q then 

lim (lMt),^(t)) = (£i,6). 
If r = and p < 0, £/ien /or any x° positive x\, 



lim (Mt),Mt)) 

t— >+oo 



(0,A//i2 



For further references, we denote by ip°° the unique point to which the system 
converges in each case. We denote by the measurable function such that ty(x°, t) 
is the value of the solution of (S r (x°) ) at time t. 
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Proof. We denote by fa and fa the functions such that (S r (x )) is written 

(2) = Zi^iW, ^(*)) and #(i) - f 2 (Mt), M*))- 

Since /i and fa are locally Lipschitz, there exists a local solution for any starting 
point x° belonging to (R + x R + )\{(0, 0)}. Moreover, for any (xi, x 2 ) <E (R+ x 
R+)\{(0,0)}, 

r - fiiXi < fi(xi,x 2 ) <r + Xpi + apx x and A(l - pi) - [i 2 x 2 < fa{xi,x 2 ) < A. 

By standard theorems about comparison of solutions of differential equations, one 
can then show that every local solution ip can be extended to R and that for any 
t e R, ip(t) = (ipi(t), ip 2 (t)) belongs to (R+ x R+)\{(0,0)}. Furthermore, with 
direct calculations, we have 

(3) + ^ t)] = r + x - kM*) - MM*)- 

For e > 0, consider 

A± = {(xi, x 2 ) e R+ x R + , < ±(r + A - ^xi - /j, 2 x 2 ) < e}, 
B\ = {(xi, x 2 ) e R+ x R + , r + A - /iiXi - [i, 2 x 2 > e}, 
B e _ = {(xi, x 2 ) E R+ x R + , r + A - fixxi - [i 2 x 2 < -s}, 

and 

A° = {(xi, x 2 ) e R+ x R + , r + A - fiiXi - [i 2 x 2 = 0}. 

According to (3), on B e + , the derivative of Vi + ^2 = \\(ipi, "02)11 is greater than 
e, hence for a starting point in A £ + , the trajectory has an L 1 increasing norm. 
Reasoning along the same lines on B e _, we see that for any 77 > 0, for any starting 
point outside ^4°, the trajectory of the differential system enters, in a finite time, 
one of the set A\ or A v _ . Moreover, upon this time, the orbit stays in the compact 
A\ U A v _ forever. It follows that (see for instance [13]) 

lim dist((V>i(t), Mt)), A ) = 0. 

This implies that any invariant set M must be included in A . We then seek 
for a maximal invariant set. It is given by the intersection of the sets Zj = 
{(xi, x 2 ), fi(xi, x 2 ) = 0}, i = 1, 2. We then remark that this system of equation 
is equivalent to the system fa + fa = and fa = 0. It turns out that 

{fa + fa){xi, x 2 ) = r + A - - ii 2 x 2 = 0. 

The equation fa{xi,x 2 ) = yields to 

li 2 x\-\x 2 

X\ = -T7Z ; -, : ^— = h(x 2 ). 

A(l - pi) - {ap + fi 2 )x 2 

The variations of h shows that ft, is a strictly decreasing diffcomorphism from / = 
[A(l — pi)/(ap + /U2), A//X2] onto R + . Hence its reciprocal function is a decreasing 
diffcomorphism from R + onto /. 

Assume first that r > 0. Then (A + r)/^i 2 > A/^2 and there exists one and only 
one equilibrium point whose coordinates (£1,^2) are thus given by the solution of 
(1) - see Figure 2 for an illustration. 

Consider the two distinct situations where \i\ = /i 2 and \x\ 7^ [i 2 . If [i\ 7^ ji 2 , then 
for any starting point (x\, x 2 ) 7^ (£1, £ 2 ) belonging to A , ^[(xi, x 2 )+ip' 2 (xi, x 2 ) = 
but (iiip[(xi, x 2 )+ [i 2 ip' 2 (xi, x 2 ) 7^ 0. Hence for t sufficiently close to 0, il>(t) does 
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Figure 2. Determination of the fixed point. 



not belong to A and then (x\, x 2 ) does not belong to M. Thus, M — {(£1, £2)} 
and according to the Poincare-Bendixson theorem (see [13] for example), 

(4) lim (lM*),^(t)) = (£i,&). 

t— > + oo 

If p\ = p 2 then ipi + ip 2 is solution of the differential equation 

v'(t) = r + A - p\v(t), u(0) = xi + x 2 . 
By direct integration, this yields to 

(Vi + lM(t) = (xi + x 2 )e-^ + —(1 - e-^). 

Mi 

This entails that A is invariant. Since A is compact, there exists a minimum 
invariant set, say M. According to the Poincare-Bendixson theorem, M is either a 
periodic orbit or a critical point. Since tpi + tp 2 is not periodic, M is also reduced 
to (£1, £2) and we have (4). 

For r = 0, the point (0, A//Z2) is a fixed point. Due to the concavity of h^ 1 , 
the sets ^4° and {/2 = 0} have at most one point of intersection with positive 
abscissa. The existence of it depends on the slope of h^ 1 at the origin. By direct 
computations, we find that 

(O'(o) = -(—+?/)■ 

M2 

Hence there exists another equilibrium point if and only if (/i _1 )'(0) > /ii//i2, i.e., 
p = ap + p 2 pi — p\ > 0. We still denote by (£1, £2) the unique solution of (1) with 
a strictly positive first coordinate. Note first that if x\ — then i\)\(t) = for any t 
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thus the vertical axis is an invariant set. Moreover, for x± — 0, a direct integration 
of (S r (x )) shows that 

lim (Mt),Mt)) = (0,A/ M2 ). 

t— > + oo 

We hereafter assume that x\ ^ 0. If ap + fi 2 pi — £*i < 0, the same reasoning as 
above shows that 

lim (Mt),Mt)) = (0, A,/ M 2). 

t— > + oo 

Assume now that ap + \i 2 pi ~ £*i > 0. At (0, X/fi 2 ), the linearization of (S r (a; )) 
gives a matrix whose determinant is given by 

d = -p^2- 

Then, according to the hypothesis, d < thus (0, X//j, 2 ) is a saddle point and 
cannot be an attractor. Reasoning as above again yields to the conclusion that 
every orbit converges to (£1, £ 2 ) f° r an y ( x i, x 2 ) such that X\ 7^ 0. □ 



5. Mean field approximation 

We now consider a sequence (X N (t) = (X^(t) 7 X 2 (t)), t > 0) of Markov pro- 
cesses with the same transitions as above but with different rates given by (with 
self evident notations): 



q?{ni, n 2 ) = r N + X N p! 



m + n 2 



q 2 (ni, n 2 ) = /Ui m 



(ni, n 2 ) = apn 2 



m + n 2 



q 4 (m, n 2 ) = X N (1 - pi 



ni + n 2 



q$ (m, n 2 ) = [i 2 n 2 . 



The main result of this Section is the following mean field approximation of the 
system X N . 



Theorem 5.1. Assume that 
T 



E 



X N (0)-x° 



> 0, — r N > r>0, — X N > X. 



Let ip(x a , .) = (ipi(x°, .), ip 2 (x° , .)) be the solution of the differential system (S r (x°) ). 
Then, for any T > 0, 



E 



sup 

t<T 



-X N (t)-^(x°,t) 



> 0. 



Before turning into the proof of Theorem 5.1, let us give the martingale problem 
satisfied by the process X N . 

Theorem 5.2. For any N > 0, the process X N is a vector-valued semi-martingale 
with decomposition: 



X»(t) =X»(0)+ f {q? + q?-q?){X ff {8))d8 + M*' 
Jo 



X 2 N (t) = X 2 (0) + / («£ - <£ q£)(X« (s)) ds + M 2 (t), 



-N / 



N N „N\i v N, 



(*) 



rN, 
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where M N = (M™ , M 2 ) is a local martingale vanishing at zero with square bracket 
given by: 

/ r* r* \ 

/ / / N , N . N\rirNi w , / N/^rN, u ■ \ 



' M 



N • 



I (^+^+^)(X N ( S )) ds - / q?(X N (s)) ds 

Jo Jo 

- f q?(X N (s)) ds [\q» + q» + q?)(X N (s)) ds 

Jo Jo 



Proof. Using the martingale problem associated with the Markov process X N , we 
get that, for t > 0, 



X N (t)=X"(0) 



N i 



( [\ q ? + q ?- q ?)(X N (s))d S 

[ ( q »- q »- q »)(X N (s))ds 
Jo 



M N (t), 



where M N = (M^ , M 2 ) is a 2-dimensional local martingale vanishing at zero. 

Let us now compute its square bracket. First of all, we consider (M^ , M 2 ). By 
integration by parts, we get that, for t > 0, 



X»{t)X»{t) = X»(0)X»(0) + f X?(s_)dX?(s) 

+ [ X»(8-)*X?{8) + [X?,X»\ t , 

J(o,t] 

where [X^ , X 2 ] denotes the mutual variation of X^ and X 2 ■ Hence 

X»(t)X»(t) = X»(0)X?(0) + fx^sM ~ 13 - a.?)(X N (s))ds 

Jo 

+ /VoO^f + q?-q?)(X N (s))ds 
Jo 

+ [X?,X?] t 

+ local martingale. 

Now, writing the martingale problem associated with the process X^ X 2 , we have 

X»(t)X?(t) = X»(0)X?(0) + fxFisM q ?)(X N (s))ds 

Jo 



+ tx?(s)( q ? - q ?)(X»(s))ds 
Jo 

+ f\x 2 N (s)-X^(s)-l)q^(X N (s))ds 
Jo 

+ local martingale. 



We conclude that 



f q${X N {s))ds. 
Jo 



Similar arguments show that 



[\q?+q?+q?)(X N (s))ds and (X 2 N ) t = [\q? +<£+<£) 
Jo Jo 



)(X N (s))ds 



which ends the proof. 



□ 
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Proof of Theorem 5.1. According to Theorem 4.1, for any x° £ R + x R + \ {(0, 0)} 
inf s6 R + ||-0(x°,s)|| > 0. Then, the theorem 5.1 is a consequence of the following 
Lemma. □ 

Lemma 5.1. There exists a constant C depending only on r, A, Pj, fii, ijl 2 and ap 
such that for any x° £ R+ x R + \ {(0, 0)}, any (X m (0))mgn sequence of random 
variables taking its values in R x R \ {(0, 0)}, for any N £ N* , and for any T > 



E 



sup 

t<T 



N 



X N (t)-^(x°,t) 



a(X (0), M £ N) 



< 



N 



— T + T 2 — \\X 



N 



N 



x cxp T / (1 + 



1 



-) 2 ds 



Proof of Lemma 5.1. Let us fix T > 0. Using Theorem 5.2, we have 

l^(t) = 1^(0) + J* + ^ - 1?)(X N (s)) ds + ±M»(t), 

±X?{t) = ±X?(0) + J* - ^ - gf)(^( S )) d, + ±M?(t). 



Moreover, 



ipi(t)= (Qi + 93 - 92) (V>(s))ds, 
Jo 

V>2(i) = / (94 - 93 - 
Jo 



g 5 ) (V'(s))ds. 



Note that for x = (x-^x 2 ) and y = (2/1,2/2) in R+ x R+ \ {(0,0)}, then 
x\ 



yi 



xi +x 2 2/1 + 2/2 



< 



xi - 2/1 



yi + 2/2 

x\ - yi 



+ 



X\ 



We also have 



< 2 



X\X2 



2/1 + 2/2 

"x-y\ 



xi + x 2 yi + 2/2 

x\ yi - X! + y 2 - x 2 



xi + x 2 



vi + yi 



1 2/11 



2/12/2 



< 2||x - 2/| 



£1 +2:2 2/1 + 2/2 

From now on, we use C for positive constants which depend only on r, A, pj, Hi, 
Ijl 2 and ap, and which may vary from line to line. For < t < T, 



lx N (t)-^(x°,t) 



(5) 
<C 



N 



X N (0)-^(x°,0) 



2 

+ T 2 


r N 


2 


VT 2 


. Aw 











+ T 



H(x°,s)\\ 



lx N (s)-^(x°,s) 



ds +J^\\ M ^1 
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Using Burkholder-Davis-Gundy inequality, we get that 

E[ sup \\M N (t)\\ 2 \a(X M (0), MeN)] < CE[| (( M w )) T | \a{X M (0), M £ N)]. 
te[o,T] 

As a consequence of Lemma 3.1 we get that for i £ {1, 2, 3, 4, 5}, 



(6) 
and 



E 



supq?(X s ) 

t<T 



<C(\\X*(0)\\ + (r N + \ N )T), 



E[\((M N )) T \\a(X M (0), M e-N)}<CT(V[\\X N (0)\\]+(r N + \ N )T). 
Hence, using Gronwall's lemma, (5) implies that 
1 



E 



sup 

t<T 



N 



X N (t)-^(x°,t) 



a(X M (0), MeN) 



< 



X N (0)-x° 



+ Ut + t 2 ±\\x n (q)\\ 



x exp 



(i 



i 



□ 



6. Stationary regime 

We have proved so far that the process N~ 1 X N converges, as N goes to infinity, 
to a deterministic R 2 -valucd function. This function converges, as t goes to infinity, 
to a fixed point tp°°. On the other hand, for each N, the Markov process X N is 
ergodic thus has a limiting distribution as t goes to infinity. This raises the natural 
question to know whether this limiting distribution converges to the Dirac mass at 
ip 00 when N goes to infinity. Let us denote by Py».„ the distribution of the process 
Y N = N~ X X N under initial distribution v. We denote by P$, v the distribution 
of the process whose initial state is chosen according to v and whose deterministic 
evolution is then given by the differential system (S r (x°j). According to Theorem 
3.1, we know that X N has a stationary probability whose value is irrespective of 
the initial distribution of X N . We denote by Y N (oo) a random variable whose 
distribution is the stationary measure of Y N . We already know that 



N- 



Y N (t),6 x o 



t—>oo 



y«(oo) 



V->oo 



U P 



t— ¥00 



The question is then to prove that this is a commutative diagram, i.e., that Y N (oo) 
converges in distribution to the Dirac measure at the equilibrium point of the 
system (S r (x°j). We borrow the proof from [15] and [7] but we need to take into 
consideration the special role of the point (0, 0) which is a singular point for some 
of the qj. 

Definition 1. We say that a probability measure v on R + x R + \ {(0, 0)} belongs 
to T when u({0, 0}) = 0. 

We will show that 1) for any sequence of initial distribution v N converging 
weakly to v with v e Vo then P y n „n converges weakly to P ,/,,„, 2) that for 
any probability measure v £ Vo, P^(t),v converges weakly to S^oo, 3) that the 
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sequence (Y N (oo), N > 1) is tight and 4) that any possible accumulation point of 
(Y N (oo), N > 1) belongs to To- 

The proof is then short and elegant: since (Y N (oo), n > 1) is tight, it is sufficient 
to prove that there is a unique possible limit to any convergent sub-sequence of 
(Y N (oo)). We still denote by Y N (oo) such a converging sub-sequence (as N goes 
to infinity). Its limit is denoted by v, known to belong to Vq. According to Point 
1. above, P y jv P converges weakly to „. Moreover by the properties of 

' y N (oo) ^ ' 

Markov processes, P y n ,p yN ; is the distribution of a stationary process, hence ip 
is also a stationary process when started from v. This means that the distribution 
of ip(t) is v for any t. Then, by Point 2. above, v = S^oo. We have thus proved 
that any convergent sub-sequence of Y N (oo) converges to 5^oo, hence the result. 
We now turn to the proof of the three necessary lemmas. 

Theorem 6.1. For any sequence of initial distribution v N converging weakly to 
v e Vq, then Pyw „w converges weakly to P^ tl/ . 

Proof. We will proceed in two steps: First prove the tightness in D([0, T], R 2 ) and 
then identify the limit. Actually, we will prove the slightly stronger result that 
Pyjv „N is tight and that the limiting process is continuous. According to [2], we 
need to show that for each positive e and ry, there exists S > and n such that for 
any N >hq, 

( \ 



sup | 

t v— u\<5 
\ v,u<T 



Y"{v)-Y"(u)\\ >e 



We denote by 



From Theorem 5.2, we know that 



Y»(v) = A?(v) + -M l N (v), i = l,2. 



Hence, for any positive a, 



(7) P sup \\Y»{v)-Y»(u)\\>e\ < P(||F iv (0)|| > a) 

\v-u\<8 
v,u<T 

+ P( sup \\A N (v) - A N (u)\\ >e/2; \\Y N (0)\\ < a) 

\v—u\<5 
v,u<T 

+ P( sup l||M»-M»||>e/2; 11^(0)11 < a). 

\v-u\<5 JV 
v,u<T 



Eqn. (6) implies that 



E 



sup^||M"( S ) 
-l n/rN 



Y N (0)\\ < a 



< 



C(a + 1) 
N ' 



This means that (N~ 1 M N , N > 1) converges to in L 2 (fl; D([0, T], R 2 ), P.|||rw(o)||<o) 
Hence it converges in distribution in D([0, T], R 2 ) and thus it is tight. This means 
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that the last summand of (7) can be made as small as needed for large N. Further- 
more, 



i=l 

fr N + \ N t C 



<2\v-u 
It follows from Lemma 3.1 that 



E 



sup \\A N {v) - A" {u)\\ 

\v— u\ <<5, 
v,u<T 



V N 



Y»(0)\\ < a 



sup \\X" ( S )|| 

s<T 



< CH r_N + *N T + a) 



< C({r + X)T + a)6. 



This means that the second summand of (7) can also be made as small as wanted. 
The hypothesis on the initial condition exactly means that this also holds for the 
first summand of (7). Thus we have proved so far that Pyn^w is tight and that 
its limit belongs to the space of continuous functions. 

We now prove that the only possible limit is P^ lV . Assume that u N tends to 
v and that V y n v n tends to some Pz,„. We suppose that the initial conditions 
X N (0) of the Markov processes are distributed as and we introduce a random 
variable x° distributed as v. Recall that Y N = N~ 1 X N . We fix M e N*, (a k = 
{a\,a%)) < k < M € R 2M + 2 and = t < h < ... < t M . We introduce 

/ r M i \ 

J2<*k,Y N (t k )> , 

-k=0 J / 



G N = E I expi 



G 



N 



E 



/ \™ x N (0) , 1\ 

I expi 2^<a k ,^{ N ,t k )> I 

V U=o ' J / 



G = E exp i 



M 



< a k ,tp(x°,t k ) > 



Lfe=o 



where Y N (-1) = 0, and Y N = N^X", with initial condition 1^(0) distributed 
as i>n and X a as v. 

Let e > 0. The sequence (z^v)jveN is tight, hence there exits a compact set 
K C R+ x R + \ {(0, 0)} such that v{K c ) + sup N v N {K°) < e. We also introduce 



Then, 



G% = E |^expi 
G K = E ^expi 
Gk = E ^exp i 



f:<a k ,Y N (t k )>] 1k(^)), 

.k=0 J / 



M 



2^ < a k ,tp( — ^,t fc ) > 



.fe=0 
' M 



TV 



i^l — i^— ) 



N 



< a k ,i>(x°,t k ) > 



.k=0 



1k(x°) . 



lim sup I G - G" I < 2e + lim sup \G% - G K \ + limsup \G K - G K \. 



N 



N N 

From Theorem 4.1, the map (x, s) h-> ip{x,s) is continuous on (R + x R + \ 
{(0,0)}) x [0,T] and inf( x iS ) e x x [o ,t] II *P( X , s)\\ > 0. Since ^-(0) takes is value in 
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the compact set K, then from Lemma 5.1, limsup^ \G^ — G^\ — 0. Since the 
sequence of measures (vn) converges weakly to v, then limsupjy \G^ — G^\ = 0. 
Hence, 

limsup|G-G w | < 2e 

N 

for all e > 0. 

That means for any to, ■ ■ ■ ,tM, 

Hence all the accumulation points are the same and the convergence of Py« v n 
follows. □ 



Theorem 6.2. For any probability measure v G Vq, P^(t),v converges weakly to 
5^oo as t — > oo. 

Proof. For any / continuous bounded on R 2 , we have 



fdP m , „ = / E \f(ij,(t)) | V(0) = x] dv(x). 
Theorem 4.1 says that for any x G R+ x R + \ {(0, 0)}, 

E \fMt)) | V(0) = x] ^ 
The result follows by dominated convergence. □ 

Theorem 6.3. The sequence (¥ N {oo), N > 1) is tight and any accumulation point 
belongs to Vq. 

We need a preliminary lemma which relies on the observation that when \i\ = 
fi2, the process X\ + Xi has the dynamics of the process counting the number 
of customers in an M/M/oo queue. Recall that = min^i, ^2) and set ( — 
(r + Xpi)/ '(J,-. For any c G R + , any x G N, define the function 

h c {t, x) = (1 + ce ^-*)* e -C C cxp( M _^ 

Note that h c is increasing with respect to x. Moreover, according to [12, Chapter 
6], 

(8) ^(t,x)+R(h c (t,.))(x)=0, 
where, for any w : N — > R, 

Rw(x) = (w(x + 1) — w(x))(r + Xpi) + (w(x — 1) — w(x))fi-X. 

Lemma 6.1. For any non negative real c, the process H c = (h c (t, Xi(t)+X 2 (t)) 7 t > 
0) is a positive supermartingale. 
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Proof. According to Dynkin formula (see [12, Proposition C.5]), for any < s < t, 
we have 



= E 



h c (t, \\X{t)\\) h c (s, \\X(a)\\) J* ^(r, ||X(r)||) 

- (r + Apj) J* (h c (r, \\X(r)\\ + 1) - h c (r, \\X(r)\\)) dr 

- f s (hc(r, \\X(r)\\ - 1) - h c (r, \\X(r)\\)) (/xi*i(r) + ^X 2 {r)) dr 
>E [m«, ||X(i)||) - M*> 11*0011) - f s ^(r, ||X(r)||) 

- (r + A PJ ) y" (fe c (r, ||X(r)|| + 1) - Mr, ll*toll)) dr 

- £(h c (r, \\X(r)\\ 1) - Mr, ||X(r) || )) ^_ (X x (r) + X 2 (r)) dr 



where the inequality follows from the monotony of h c and the definition of jU_ 
Hence we get that 



> E 



h c (t,\\X(t)\\)-h c (s,\\X(s)\\) 



dh. 



-^(r, ||X(r)||) + i2(Mr,.))(ll^(r)||)dr 

In view of Eqn. (8), we get 

0>-E[h c (t,\\X(t)\\)- h c (s,\\X(s)\\)\T s ], 

i.e., H c is a supermartingalc. 

Now, let (Y N (oo), N > 1) be a subsequence which converge to v. Since X w (oo) 
is a random variable distributed according to the stationary law of the process X N , 



E 



0. 



Qe-"(X N (oo)) 
By a direct calculation, we have 

Qe-"(x) = e~ W [(A + r)^ 1 - 1) + faxi + /x 2 z 2 )(e - 1)], 

then 

-V||V«(oo) 



(AA/+r A r)(l-e- 1 )E 



-V||V N (oo) 



V/ 



= JV(e-l)E 

Hence, (1 - e" 1 )^ + A)z/({(0, 0)}) = 0, i.e., v belongs to Po 

Proof of Theorem 6.3. Let X be real, for any positive real 9, we have 

Pfliy^t)!! > K) = P(\\X N (t)\\ > NK) < e- eNK E [exp^X^)!!)] . 
Lemma 6.1 entails that 
E [exp^UX^)!!)] < (1 + {e e - i) e -f-t)JV|l* w (o)ll exp^A^e 9 - 1)(1 - e""-*)) . 



□ 
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Hence. 



P(\\Y»(oo)\\>K)=Vm P(\\Y"(t)\\>K) 

t— ¥00 

_ t -.N\\X N (0)\\ 



f llX m exp(-6NK + N((e e - 1)(1 - e""-*)) 



< inf lim fl + (e - l)e" 

6»>0t^oo v 

= inf exp (N(-6K + C{e e - 1)) 

< exp(-^NK In K), 

for K large enough. The tightness follows. 

7. Central Limit Theorem 



□ 



It turns out that we can also evaluate the order of the approximation when we 
replace X N by ip. This is given by CLT like theorem. 

Theorem 7.1. Assume that the hypothesis of Theorem 5.1 holds. Then, for any 
T > 0, the process 

w N = Vn(y n - v) 

tends in distribution in D([0, T], R 2 ) to a centered Gaussian process with covariance 
matrix T(t) given by: 

( 



T{t) 



-up / — 

Jo i>i 



ri(t) 



-ap 



o ipi(s) + ip 2 (s) 



where 



(s) +1p2(s 

Y l{ t)=rt + [\p I - 
Jo V>i 

r 2 (t)= f\{i 

Jo 



ds 



r 2 (t) 



/ 



Vi(s) 



Hi Vl(S) + ap ; ; ; , ; ; OS 



) + H2V2(s) +ap 



M*) + M*) 

4>i(s)ip 2 (s) 



tpl(s) + tp 2 (s) 1pl(s) + l/>2(s) 

Proof. According to [5, p. 339], it suffices to prove that 



E 



BUp\W"(t)-W N (t-)\ 
t<T 



JV->+oo 



> 



and that 



'W N h^^T(t). 



Since the jumps of Y N are bounded by 1/N, those of W N are bounded by N 1//2 , 
hence the first point is proved. As to the second point, remark that 



W N )) t 



N-H(M N ' 



and then use Theorem 5.1. 



□ 



8. Numerical investigation 

Another approach to evaluate the order of approximation can be made by com- 
puter simulation. We simulated the Markov process for N — 100 and compute the 
estimate of the prevalence by a simple Monte-Carlo method on 10, 000 trajectories. 
For the parameters we chose, a = 1, /j,i = 0.1, /j,2 = 0.2, r = 1, A = 5 and pi = 0.8, 
the results are strikingly good as shown in Figure 3. Note that the choice of pa- 
rameters is here very delicate since biological parameters arc not very well known 
(i.e., Hi, H2, P, • • •) and population dependant quantities are even more obscure to 
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D Ol I1Z 0.3 O.A Z>'. 

Figure 3. Prevalence with respect to p. The solid line represents 
the value as computed by Equations 1. The dots represents the 
simulated values. The 95% confidence interval are so small, they 
can't be displayed. 



determine. We here chose parameters which seems reasonable and fit the observed 
prevalence. 

In such models, another quantity of interest is the relative importance of each 
parameters: what does affect most the prevalence ? On the deterministic system, 
this question is easily solved by computing the derivative of the prevalence with 
respect to each of the parameters. We now explain how to compute the sensitivity 
of the prevalence on the stochastic model. Say we have a function F bounded which 
depends on the sample-paths of X, we aim to compute: 



dp 



E P tn 



where we put a p under the expectation symbol to emphasize the dependence of 
the underlying probability with respect to p. Other "greeks", as these quantities 
are called in mathematical finance, can be derived analogously. We assume that 
we observe the Markov process on a time window of size T, i.e., any functional is 
implicitly assumed to belong to Ft = u{X(s), < s < T}. 



Theorem 8.1. For any bounded F, F £ Ft, we have: 



dp 



F 



[£i{(i,-i)}(A*(s))- f 

\s<T Jo 



q 3 (X(s-)) ds 



-con,, I F, ^l {(1) _ 1)} (AJf(s)) 

s<T 



where AX(s) = X{s) - X(s_). 

Proof. The proof relies on the Girsanov theorem which is more easily expressed 
in the framework of multivariate point measures. Since there are only five kind of 
jumps, we can represent the dynamics of X as a point measures on R + x {1, • • • , 5}: 



H([0, t] x {*}) = 



^2l{AX(s)=h}, 

s<t 
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where 



h = (1, 0), h = (-1, 0), ; 3 = (1, -1), h = (0, 1) and l 5 = (0, -1). 



In the reverse direction, 



X(t) = X(0) + j>([0, t]x 



It is immediate from the preceding results that v v , the P p -predictable compensator 
of is given by 

diA>(t, i) = qi (X(t_))dt. 
To compute d/dpE p [F] means to compute 



lim-(E p+£ [F]-E p [F]). 

s— >0 e 



Under P p+e , 



di/ p+e (t, i) = dv p (t, i) for i £ 3 and d^ p+£ (t, 3) = (1 + -) di/ p (t, 3). 



Let 



U(t,i)={e_ . f . = 3> 



ifi^3, 

According to the Girsanov theorem (see [3, 8]), this means that 



E P + £ [F] — E p 



E„ 



F£( [ U{s, i)(dfx(s, i)-dv p (s, i))) 
Jo 



F£(-(»([0, T]x{3})-^([0, T]x{3}))), 



where £ denotes the Doleans-Dade exponential. It is known (see [3]) that a Doleans- 
Dade exponential follows the same rule of derivation as a usual exponential, hence 
the result. □ 



With the parameters above, the simulated greek coincides pretty well with the 
sensitivity computed by differentiating the expression of the stationary prevalence 
in the deterministic system, see Figure 4. However, as usual with this method, the 
confidence interval are rather large. 
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>:• 41 or 1X3 oi as o& or 

p 

Figure 4. Prevalence greek with respect to p. Same conventions 
as above. The error bars represent the 95% confidence interval. 
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